import numpy as np
import csv
from datetime import datetime
from time import strftime


def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def build_data_list(inputCSV):
    sKey = []
    fn = inputCSV
    ra = csv.DictReader(file(fn), dialect="excel")
    
    for record in ra:
        #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
        for item in ra.fieldnames:
            temp = float(record[item])
            sKey.append(temp)
    sKey = np.array(sKey)
    sKey.shape=(-1,len(ra.fieldnames))
    return sKey

def calDistance(tList):
    tDist = []
    for item in tList:
        temp = calItemListDist(item, tList)
        tDist.append(temp)
    tDist = np.array(tDist)
    tDist.shape = (len(tList), -1)
    return tDist

def calItemListDist(tItem, tList):
    tDist = []
    for item in tList:
        temp = ((tItem[0] - item[0])**2 + (tItem[1] - item[1])**2)**0.5
        tDist.append(temp)
    return tDist

def kernel_weight(data, bandwidth = None, fixed = True, k = 2, function = 'triangular'):
    '''Spatial weights based on kernel functions
    http://pysal.org/1.1/library/weights/Distance.html#pysal.weights.Distance.Kernel
    
    Parameters
    ----------

    data        : array (n,k)
                  n observations on k characteristics used to measure
                  distances between the n objects
    bandwidth   : float or array-like (optional)
                  the bandwidth :math:`h_i` for the kernel. 
    fixed       : binary
                  If true then :math:`h_i=h \\forall i`. If false then
                  bandwidth is adaptive across observations.
    k           : int
                  the number of nearest neighbors to use for determining
                  bandwidth. For fixed bandwidth, :math:`h_i=max(dknn) \\forall i`
                  where :math:`dknn` is a vector of k-nearest neighbor
                  distances (the distance to the kth nearest neighbor for each
                  observation).  For adaptive bandwidths, :math:`h_i=dknn_i`
    function    : string {'triangular','uniform','quadratic','quartic','gaussian'}
                  kernel function defined as follows with 

                  .. math::

                      z_{i,j} = d_{i,j}/h_i

                  triangular 

                  .. math::

                      K(z) = (1 - |z|) \ if |z| \le 1

                  uniform 

                  .. math::

                      K(z) = |z| \ if |z| \le 1

                  quadratic 

                  .. math::

                      K(z) = (3/4)(1-z^2) \ if |z| \le 1

                  quartic

                  .. math::

                      K(z) = (15/16)(1-z^2)^2 \ if |z| \le 1
                 
                  gaussian

                  .. math::

                      K(z) = (2\pi)^{(-1/2)} exp(-z^2 / 2)
    '''
    if k > (len(data) - 1):
        print k, ' exits the boundary of data'
        return
    distance = calDistance(data)
    if bandwidth == None:
        if fixed:
            # if fixed == True
            tDist = np.sort(distance, axis=1)
            tMax = np.max(tDist[:,k])
            bandwidth = np.ones((len(data),1),'float') * tMax
        else:
            tDist = np.sort(distance, axis=1)
            bandwidth = tDist[:,k]

    kernel = np.zeros((len(data),len(data)),'float')

    c = np.pi*2
    c = c**(-0.5)
    i = 0
    for item in distance:
        j = 0
        for tIDist in item:
            if tIDist <= bandwidth[i]:
                z = tIDist/bandwidth[i]
                if function == 'triangular':
                    kernel[i,j] = 1 - z
                elif function == 'uniform':
                    kernel[i,j] = z
                elif function == 'quadratic':
                    kernel[i,j] = (3./4)*(1-z**2)
                elif function == 'quartic':
                    kernel[i,j] = (15./16)*(1-z**2)**2
                elif function == 'gaussian':
                    kernel[i,j] = c*np.exp(-(z**2)/2.) 
                else:
                    print 'Unsupported kernel function', function
                    return
            j += 1                
        i += 1
    return kernel

def kernel_smoother(e, b, weight):
    kRate = []
    for i in range(0, len(e)):
        temp_e = 0.0
        temp_b = 0.0
        for j in range(0, len(e)):
            temp_e += e[j] * weight[i,j]
            temp_b += b[j] * weight[i,j]
        kRate.append(temp_e/temp_b)
    return kRate

def test():
    points = [(10, 10), (20, 10), (40, 10), (15, 20), (30, 20), (30, 30)]
    e = np.array([30, 25, 25, 15, 33, 21])
    b = np.array([100, 100, 110, 90, 100, 90])
    mykw = kernel_weight(points)
    temp_rate = kernel_smoother(e, b, mykw)
    print temp_rate

#--------------------------------------------------------------------------
#MAIN
if __name__ == "__main__":
    print "begin at " + getCurTime()
    #test()
    unitCSV = 'C:/_DATA/CancerData/SatScan/mult6000/three16_format.csv'
    geoCSV = 'C:/_DATA/CancerData/SatScan/neast_id_geo.csv'
    dataMatrix = build_data_list(geoCSV)
    points = dataMatrix[:, 1:3]
    mykw = kernel_weight(points, k = 5, fixed = False, function='gaussian')
    dataMatrix = build_data_list(unitCSV)
    pop = dataMatrix[:, 1]
    kernel_rate = []
    for i in range(0,1):
        print i
        cancer = dataMatrix[:, 2+i]
        temp_rate = kernel_smoother(cancer, pop, mykw)
        kernel_rate.append(temp_rate)
    #kernel_rate = np.array(kernel_rate)
    #kernel_rate.shape = (1000, -1)
    #kernel_rate = np.transpose(kernel_rate)
    print kernel_rate

    filePath = 'C:/_DATA/CancerData/SatScan/mult6000/three16_kernel_rate.csv'
    #np.savetxt(filePath, kernel_rate, delimiter=',')
    print "end at " + getCurTime()
    print "==========================="